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Impact  of  Assimilating  Surface  Velocity  Observations  on  the  Model  Sea 
Surface  Height  Using  the  NCOM-4DVAR* 

Matthew  J.  Carrier,  Hans  E.  Ngodock,  Philip  Muscarella,  and  Scott  Smith 

Naval  Research  Laboratory,  Stennis  Space  Center,  Mississippi 

(Manuscript  received  8  September  2014,  in  final  form  8  December  2015) 

ABSTRACT 

The  assimilation  of  surface  velocity  observations  and  their  impact  on  the  model  sea  surface  height  (SSH)  is 
examined  using  an  operational  regional  ocean  model  and  its  four-dimensional  variational  data  assimilation 
(4D  VAR)  analysis  component.  In  this  work,  drifter- derived  surface  velocity  observations  are  assimilated  into 
the  Navy’s  Coastal  Ocean  Model  (NCOM)  4DVAR  in  weak-constraint  mode  for  a  Gulf  of  Mexico  (GoM) 
experiment  during  August-September  2012.  During  this  period  the  model  is  trained  by  assimilating  surface 
velocity  observations  (in  a  series  of  96-h  assimilation  windows),  which  is  followed  by  a  30-day  forecast  through 
the  month  of  October  2012.  A  free-run  model  and  a  model  that  assimilates  along-track  SSH  observations  are 
also  run  as  baseline  experiments  to  which  the  other  experiments  are  compared.  It  is  shown  here  that  the 
assimilation  of  surface  velocity  measurements  has  a  substantial  impact  on  improving  the  model  representa¬ 
tion  of  the  forecast  SSH  on  par  with  the  experiment  that  assimilates  along-track  SSH  observations  directly. 
Finally,  an  assimilation  experiment  is  done  where  both  along-track  SSH  and  velocity  observations  are  utilized 
in  an  attempt  to  determine  if  the  observation  types  are  redundant  or  complementary.  It  is  found  that  the 
combination  of  observations  provides  the  best  SSH  forecast,  in  terms  of  the  fit  to  observations,  when  com¬ 
pared  to  the  previous  experiments. 


1.  Introduction 

Accurate  representation  of  the  sea  surface  height  is  an 
important  aspect  of  any  operational  ocean  forecasting 
system.  In  the  deep  ocean,  and  in  the  absence  of  strong 
surface  wind  forcing,  the  ocean  is  in  quasigeostrophic 
balance  where  the  gradient  of  the  SSH  defines  the  sur¬ 
face  currents.  Events  such  as  surface  oil  spills,  tracking 
of  downed  aircraft,  or  searching  for  lost  ships  or  persons 
at  sea  are  directly  affected  by  surface  ocean  currents. 
Accurate  prediction  of  SSH  can  be  difficult,  however,  as 
ocean  mesoscale  eddies  (which  account  for  much  of  the 
ocean  SSH  variability)  are  nondeterministic  and  small 
errors  in  the  model  grow  over  time  rendering  lengthy 
prediction  impossible  without  regular  correction  by  the 
assimilation  of  observations  (Jacobs  et  al.  2014).  Ocean 


*  Naval  Research  Laboratory  Contribution  Number  NRL/JA/ 
7320-14-2292. 


Corresponding  author  address :  Matthew  J.  Carrier,  Naval  Re¬ 
search  Laboratory,  Bldg.  1009  Batch  Blvd.,  Stennis  Space  Center, 
MS  39529. 

E-mail:  matthew.carrier@nrlssc.navy.mil 


circulation  models,  such  as  the  Navy  Coastal  Ocean 
Model  (NCOM;  Martin  2000;  Barron  et  al.  2006),  the 
Hybrid  Coordinate  Ocean  Model  (HYCOM;  Bleck 
2002),  Princeton  Ocean  Model  (POM;  Blumberg  and 
Mellor  1987),  and  the  Regional  Ocean  Modeling  System 
(ROMS;  Shchepetkin  and  McWilliams  2003,  2005; 
Marchesiello  et  al.  2001),  among  others,  use  widely  avail¬ 
able  sea  surface  height  anomaly  (SSHA)  measure¬ 
ments  from  satellite  altimeters  to  regularly  correct  the 
model  SSH. 

There  exist  a  number  of  methods  for  assimilating 
satellite  altimeter  data  into  ocean  models.  One  general 
method  involves  adjusting  the  subsurface  thermody¬ 
namic  structure  of  the  model  (i.e.,  Haines  1991;  Cooper 
and  Haines  1996;  Fox  et  al.  2002)  in  order  to  generate 
adjustments  to  the  SSH.  These  methods  are  mainly 
employed  with  less  sophisticated  data  assimilation  (DA) 
methods,  such  as  optimal  interpolation  (OI)  or  three- 
dimensional  variational  data  assimilation  (3D VAR), 
when  a  multivariate  background  error  covariance  is  not 
easily  attainable.  More  advanced  methods  that  are  ca¬ 
pable  of  providing  a  multivariate  background  error  co- 
variance,  either  directly  or  through  the  action  of  an 
adjoint  (AD)  and  tangent  linear  (TL)  model,  are 
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capable  of  assimilating  altimeter  data  without  resorting 
to  empirical  derivation  of  subsurface  structures.  In  these 
cases,  the  altimeter  data  are  either  assimilated  as  an 
interpolated  map  product,  such  as  in  Moore  et  al.  (2011), 
Zavala-Garay  et  al.  (2012),  or  Xu  et  al.  (2013),  among 
others,  or  as  a  direct  assimilation  of  along-track  obser¬ 
vations,  as  in  Hoteit  et  al.  (2013).  Regardless  of  whether 
one  uses  along-track  or  interpolated  data,  it  remains  that 
one  must  convert  the  anomaly  measurement  to  the  form 
of  the  ocean  model  in  order  to  assimilate  the  data.  Doing 
so  requires  knowledge  of  the  geoid,  which  can  only  be 
estimated,  to  convert  the  anomaly  information  into  dy¬ 
namic  height.  In  practice,  a  long-term  mean  from  ob¬ 
servations  or  a  numerical  model  is  used  as  a  proxy  for 
the  geoid  (Chassignet  et  al.  2007). 

Ocean  velocity  observations  have  the  potential  to 
provide  information  regarding  the  shape  and  gradient  of 
the  ocean  SSH,  without  the  need  to  derive  an  estimate  of 
the  geoid.  Observing  systems,  such  as  Acoustic  Doppler 
Current  Profilers  (ADCP),  high-frequency  (HF)  radar, 
and  surface  ocean  drifters,  such  as  NO  A  A  drifting 
buoys,  can  provide  information  on  ocean  currents.  Let 
us  consider  the  case  of  surface  drifters.  It  is  well  known 
that  surface  currents  are  generally  in  quasigeostrophic 
balance  with  the  dynamic  ocean  height  field,  except  in 
cases  of  strong  surface  wind  forcing,  or  near  the  coast¬ 
lines  where  river  inflow  and  offshore  winds  play  a  role. 
Therefore,  in  the  absence  of  strong  surface  winds  and  for 
locations  far  from  shore,  a  passive  surface  drifter  will 
follow  the  contours  of  the  SSH  and  the  speed  of  the 
drifter  provides  information  regarding  the  gradient  of 
this  field.  The  question  becomes,  then,  can  surface  ve¬ 
locity  observations  help  to  constrain  the  model  SSH  field 
in  the  absence  of  along-track  SSH  observations?  Fur¬ 
ther,  can  along-track  SSH  observations  be  combined 
with  in  situ  surface  velocity  data  to  further  improve  the 
model  forecast  of  SSH? 

A  recent  experiment,  dubbed  the  Grand  Lagrangian 
Deployment  (GLAD),  was  conducted  in  the  Gulf  of 
Mexico  (GoM)  in  the  latter  half  of  2012  by  the  Con¬ 
sortium  for  Advanced  Research  on  Transport  of  Hy¬ 
drocarbon  in  the  Environment  (CARTHE).  Roughly 
300  Coastal  Ocean  Dynamics  Experiment  (CODE)- 
type  surface  drifters  (drogued  to  1  m)  were  released  in 
July  2012,  and  their  positions  were  tracked  through  the 
GoM  for  the  remainder  of  the  year  (Poje  et  al.  2014). 
These  drifters  reported  their  positions  every  5  minutes, 
making  the  estimation  of  accurate  Eulerian  velocities 
possible.  An  earlier  effort  by  Carrier  et  al.  (2014)  uti¬ 
lized  the  Navy’s  Coastal  Ocean  Model  (NCOM) 
4DVAR  analysis  system  to  assimilate  Eulerian  velocity 
estimates  derived  from  the  GLAD  drifters  with  NCOM 
forecasts  of  the  GoM.  The  results  showed  substantial 


improvement  in  the  model  forecast  of  surface  currents  in 
the  vicinity  of  the  GLAD  drifters,  with  skill  out  to  96  h. 
In  doing  so,  Carrier  et  al.  (2014)  demonstrated  the  fea¬ 
sibility  of  surface  drifter-derived  velocity  assimilation 
and  prediction.  In  this  current  work,  we  perform  addi¬ 
tional  experiments  and  expand  the  analysis  begun  in 
Carrier  et  al.  (2014)  to  quantitatively  assess  the  impact 
of  velocity  assimilation  on  the  model  SSH.  To  that  end, 
four  experiments  are  presented  in  this  work:  1)  an  as¬ 
similation  of  GLAD-derived  surface  velocity  observa¬ 
tions,  2)  an  assimilation  of  along-track  SSH  observations, 
3)  a  combined  assimilation  of  GLAD-derived  velocities 
and  along-track  SSH,  and  4)  a  nonassimilative  free  run. 
The  results  from  these  experiments  are  examined  and 
compared  to  observations. 

In  the  next  section  the  model,  assimilation  system,  and 
the  observations  are  introduced;  section  3  presents  an 
overview  of  the  GoM  circulation  during  the  experiment 
time  frame;  and  the  experiment  design,  results,  and 
conclusions  are  given  in  section  4,  followed  by  a  sum¬ 
mary  in  section  5. 

2.  Dynamical  model,  analysis  system,  and 
observations 

a.  Dynamical  ocean  model 

The  dynamical  model  used  for  this  work  is  NCOM, 
which  is  a  primitive  equation  ocean  model  with  a  free 
surface  and  a  generalized  vertical  coordinate  that  can  be 
configured  with  terrain-following  free  sigma  or  fixed 
sigma,  or  constant  z-level  surfaces  in  a  number  of  com¬ 
binations  (Barron  et  al.  2006).  The  model  employs  the 
Mellor-Yamada  level-2.5  turbulence  closure  parame¬ 
terization  (Mellor  and  Yamada  1982)  for  vertical  diffu¬ 
sion  and  the  Smagorinsky  scheme  (Smagorinsky  1963) 
for  horizontal  diffusion. 

The  model  domain  for  all  experiments  extends  from 
79°-98°W  to  18°-31°N  using  a  spherical  coordinate 
projection  at  a  horizontal  resolution  of  6  km.  The  model 
has  50  layers  in  the  vertical  extending  down  to  a  maxi¬ 
mum  of  5500  m.  Lateral  boundary  conditions  for  each 
experiment  are  provided  by  the  global  NCOM  at  Vs0 
resolution;  global  NCOM  does  not  include  the  effect  of 
tides.  For  this  study  the  24-h  forecast  from  global 
NCOM,  valid  on  1  August  2012,  is  used  to  provide  the 
initial  condition;  this  is  done  for  the  very  reason  that  it  is 
not  as  good  as  the  analysis  field,  but  is  very  close  to  it,  so 
that  the  improvements  gained  from  the  data  assimilation 
will  be  more  apparent.  It  should  be  noted  that  the  global 
NCOM  model  has  the  same  vertical  structure  as  the 
regional  model  used  in  this  study  (only  the  resolution  in 
the  horizontal  is  different).  The  global  NCOM  was,  as  of 
August  2012,  the  operational  global  forecasting  model 
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used  by  the  U.S.  Navy  (Barron  et  al.  2006).  It  was 
an  assimilative  model,  using  the  multivariate  optimal 
interpolation  (MVOI)  Navy  Coupled  Ocean  Data  As¬ 
similation  (NCODA)  system  (Cummings  2005).  Anal¬ 
ysis  updates  were  generated  every  24  h  using  observations 
in  a  ±  12-h  window  around  the  analysis  time.  Surface  at¬ 
mospheric  forcing,  such  as  wind  stress,  atmospheric 
pressure,  and  surface  heat  flux  is  provided  by  the  0.5° 
NOGAPS  model  every  3h  (Rosmond  et  al.  2002);  river 
forcing  is  provided  via  an  internal  NRL  river  product  that 
includes  monthly  mean  river  data  for  each  major  river 
across  the  globe.  There  is  no  tidal  forcing  added  to  the 
boundary  conditions  for  these  experiments. 

b.  Ocean  analysis  system 

The  DA  system  selected  for  this  work  is  the  NCOM- 
4DVAR  developed  using  the  dynamical  core  of  NCOM. 
The  NCOM-4DVAR  is  a  variational  assimilation  sys¬ 
tem  based  on  the  indirect  representer  method  as  de¬ 
scribed  by  Bennett  (1992,  2002)  and  Chua  and  Bennett 
(2001).  This  system  has  been  described  in  detail  by 
Ngodock  and  Carrier  (2013)  and  Ngodock  and  Carrier 
(2014),  therefore,  only  a  brief  description  is  provided 
here.  The  representer  method  aims  to  find  an  optimal 
analysis  solution  as  the  linear  combination  of  a  first 
guess  (i.e.,  prior  model  solution)  and  a  finite  number  of 
representer  functions,  one  per  datum: 

M 

u(x,  t)  =  uF(x,  t)  +  X  Pmrm(x,  t),  (1) 

m= 1 

where  u(x ,  t)  is  the  optimal  analysis  solution,  uF(x ,  t)  is 
the  prior  forecast,  rm(x,  t )  is  the  representer  function  for 
the  rath  observation,  and  is  the  rath  representer  co¬ 
efficient.  The  representer  coefficients  can  be  found  by 
solving  the  linear  system: 

(R  +  0)/3  =  y  —  Hx^ ,  (2) 


until  convergence.  Once  found,  /3m  is  acted  upon  in  (1), 
involving  one  final  sweep  through  the  adjoint  and  TL 
models  to  find  the  analysis  increment. 

The  background  and  model  error  covariance  in 
NCOM-4DVAR  follow  the  work  of  Weaver  and 
Courtier  (2001)  and  Carrier  and  Ngodock  (2010),  where 
the  error  correlation  portion  of  the  covariance,  for  both 
the  model  and  the  background  errors,  are  not  directly 
calculated  and  stored  in  NCOM-4DVAR;  rather,  the 
effect  of  the  correlation  matrix  acting  on  an  input  vector 
is  modeled  by  the  solution  of  a  diffusion  equation.  For  a 
full  explanation  of  this  method,  we  refer  the  reader  to 
Weaver  and  Courtier  (2001)  or  Yaremchuk  et  al.  (2013); 
for  a  description  of  the  implementation  of  this  method  in 
NCOM-4DVAR,  we  refer  the  reader  to  Carrier  and 
Ngodock  (2010)  or  Ngodock  (2005). 

The  multivariate  and  anisotropic  characteristics  of 
the  4DVAR  error  covariance  are  achieved  by  the  lin¬ 
ear  dynamical  balance  relationships  that  are  part  of  the 
AD  and  TL  models.  Methods  based  on  3DVAR  at¬ 
tempt  to  model  the  multivariate  error  covariance  by 
way  of  empirical  balance  operators  (Weaver  et  al. 
2005),  whereas  ensemble  methods  rely  on  the  ensemble 
covariance  itself  to  provide  these  connections  (Ngodock 
et  al.  2006).  In  the  case  of  how  velocity  observations 
can  impact  the  model  SSH,  the  NCOM-4DVAR  has 
two  primary  mechanisms:  geostrophic  balance  and 
continuity. 

For  the  geostrophic  balance  we  can  examine  the 
NCOM  momentum  equations,  which  are  expressed  as 


du  _  ,  x  _  1  dp  d  ( du 

dt  y  j  ^  j  podx  u  dz\  Mdz 

and 


(3) 


where  O  is  the  observation  error  covariance,  y  is  the 
observation  vector,  H  is  the  linear  observation  operator 
that  maps  the  model  fields  to  the  observation  locations, 
is  the  model  vector,  and  R  is  the  representer  matrix 
and  is  equivalent  to  HMBMTHT  (M  is  the  TL  model;  MT 
is  the  adjoint  model;  B  is  the  initial  condition  or  model 
error  covariance,  depending  on  what  portion  of  the  y  - 
Hx^  vector  it  is  applied  to;  and  the  superscript  T  denotes 
the  linear  transposition).  Since  the  matrix  R  +  O  is 
symmetric  and  positive  definite,  (2)  can  be  solved  for 
P  iteratively  using  a  linear  solver,  such  as  the  conjugate 
gradient  method.  From  (2)  it  is  clear  that  the  f3m  for  each 
representer  can  be  found  by  integrating  the  adjoint  and 
TL  models  over  some  number  of  minimization  steps 


where  u  and  v  are  the  velocity  components,  Q  is  the 
volume  flux  source  term,  /is  the  Coriolis  parameter,  p  is 
the  pressure,  pQ  is  the  reference  water  density,  F  is  the 
horizontal  mixing  term  for  momentum,  and  Km  is  the 
vertical  eddy  coefficient  for  momentum.  Let  us  consider 
just  the  horizontal  pressure  gradient  component  of  (3) 
and  take  its  first-order  derivative: 


d8u 

1 

d8p  + 

~dT~ 

Po 

dx 

d8v 

1 

d_8p_ 

Po 

dy 

with  its  adjoint  as 
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Fig.  1.  (a)  NCOM  Gulf  of  Mexico  domain  with  location  of  Dirac  Delta  impulse  indicated  by 
red  star  and  (b)  SSH  cross  correlation  (color)  with  velocity  perturbations  (arrows)  overlaid  at 
end  of  96-h  TL  model  integration. 


dA  1  /3A  3A  \ 

— p-  =  —  — ±  +  — G 

dt  po  V  dx  dy) 


(5) 


It  is  clear  from  (5)  that  a  forcing  on  the  adjoint  velocity 
field  (A„  and  A,,)  will  result  in  an  update  for  the  adjoint 
pressure  variable  \p.  The  horizontal  gradient  of  the 
pressure  field  is  defined  in  NCOM  as 


g«+s  f°EJz  and 

po  dx  po  dx  dx  po  Jz  dx 

1  Sp_  1  dp(Q  ,  d£  ,  g 
PQ  dy  Po  dy  g  dy  po  )z  dy 


(6) 


where  f  is  the  surface  elevation  and  g  is  gravity.  Exam¬ 
ining  the  components  involved  with  the  horizontal 
pressure  gradient  due  to  differences  in  the  surface  ele¬ 
vation  and  taking  its  first-order  derivative  gives 


1  d8p 
p  dx 

r  o 

1  d8p 

Po  dy 


d8£ 

...  +  ... 

5  dx 


and 


-g—  + 

dy 


(7) 


and  its  adjoint  is  given  as 


dx 

d\g 

dy 


-P08 


-Po8 


d\ 

_ p_ 

dx 

d\ 

P 

dy 


and 


(8) 


Therefore,  from  (5)  and  (8),  a  forcing  on  adjoint  velocity 
will  propagate  to  the  adjoint  of  the  pressure  gradient,  which 
in  turn  will  force  the  adjoint  of  the  surface  elevation  (i.e., 
SSH)  via  the  adjoint  of  the  geostrophic  balance  relation¬ 
ship.  In  practice,  the  adjoint  solution  is  used  to  initialize  the 
TL  model,  and  therefore  the  TL  SSH  field  receives  this 
perturbation  and  propagates  it  forward  in  time. 

The  TL  model  provides  its  own  mechanism  for  ve¬ 
locity  observation  information  to  influence  the  SSH  field 


during  the  forward  integration  of  the  NCOM-4DVAR. 
In  NCOM  the  depth-integrated  momentum  and  conti¬ 
nuity  equations  are  needed  to  calculate  the  surface  el¬ 
evation.  The  depth-integrated  continuity  equation  is 
given  as 


d£  =  _d(Du)  _  a(Du)  +  Dq 
dt  dx  dy 

where  D  is  the  total  depth,  Q  is  the  depth-integrated 
volume  flux  source  term,  and  u  and  v  are  depth- 
integrated  baroclinic  velocity  terms.  The  linearization 
of  (9)  is  straightforward  and  indicates  that  the  TL  ve¬ 
locity  components  will  induce  a  perturbation  on  the  TL 
surface  elevation.  It  is  through  these  dynamical  balance 
relationships,  based  on  the  nonlinear  NCOM,  that  the 
NCOM  TL  and  AD  models  are  able  to  extract  accurate 
sea  surface  height  adjustments  from  surface  velocity 
observations. 

One  can  examine  this  feature  in  isolation  by  per¬ 
forming  one  run  of  the  AD  and  TL  models  using  a  pair 
of  impulses  on  one  set  of  the  adjoint  velocity  compo¬ 
nents.  These  impulses  are  defined  as  Dirac  delta  func¬ 
tions  centered  at  the  observation  locations,  or  a  single 
point  in  the  model  domain.  The  AD  model  is  then  in¬ 
tegrated  backward  to  the  initial  time,  the  result  is  con¬ 
volved  by  the  appropriate  covariance  functions  and  is 
used  to  force  the  TL  model  forward  to  the  end  of  the 
integration  time.  Ligure  la  shows  a  map  of  the  NCOM 
GoM  domain.  The  red  star  indicates  the  location  of  a 
Dirac  impulse  forcing  to  the  adjoint  model;  in  this  case, 
Dirac  impulses  are  placed  on  the  u-  and  u-velocity 
components.  Figure  lb  shows  the  resulting  cross  corre¬ 
lation  of  velocity  to  SSH  (with  velocity  perturbations 
overlaid)  at  the  end  of  the  integration  window,  which  in 
this  case  is  96  h.  The  cross  correlation  to  SSH  shown  here 
is  strictly  due  to  the  propagation  of  the  velocity  in¬ 
formation  through  the  linearized  model  dynamics;  there 
is  no  external  forcing  (surface,  boundary)  added  to  the 
TL  and  AD  models.  Since  the  Dirac  impulses  placed  on 
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the  velocity  components  were  both  positive  unit  values, 
the  resulting  velocity  field  is  oriented  toward  the 
northeast  with  a  perturbation  on  SSH  that  forms  a  di¬ 
pole  structure,  with  a  positive  (negative)  height  anomaly 
to  the  southeast  (northwest)  of  the  initial  impulse  loca¬ 
tion.  The  asymmetric  appearance  of  the  structure  is  due 
to  the  influence  of  the  background  state  (around  which 
the  TL  model  is  linearized).  A  similar  pattern  can  be 
generated  by  employing  a  Dirac  impulse  on  the  SSH 
variable  (not  shown). 

c.  GLAD  and  SSH  observations 

The  observations  used  in  this  study  are  the  GLAD- 
derived  surface  velocity  observations  as  well  as  altimeter 
derived  estimates  of  the  sea  surface  height.  The  deri¬ 
vation  of  surface  velocities  from  the  GLAD  data  in¬ 
volves  the  calculation  of  the  position  difference  of  the 
drifter  over  some  time  scale.  In  practice,  however,  the 
calculation  is  more  complicated  due  to  noisy  position 
values  and  nonphysical  accelerations.  Because  of  this,  a 
multistep  process  is  used  to  derive  Eulerian  velocities 
from  the  GLAD  drifter  position  data.  The  details  of  the 
processing  can  be  found  in  Carrier  et  al.  (2014)  and 
Muscarella  et  al.  (2015).  The  result  of  this  processing  is  a 
dataset  of  surface  velocity  observations  available  in 
15-min  intervals  for  each  drifter.  The  observations  are 
given  an  error  standard  deviation  value  of  0.02 ms-1 
(Carrier  et  al.  2014). 

Altimeter  data  are  obtained  from  an  array  of  polar- 
orbiting  satellites.  These  data  are  processed  through 
the  Altimeter  Processing  System  (ALPS;  Jacobs  et  al. 
2002),  which  is  available  from  the  Altimetry  Data  Fu¬ 
sion  Center  (ADFC)  at  the  Naval  Oceanographic  Of¬ 
fice  (NAVOCEANO).  These  data  are  stored  as  SSHA 
values  and  must  be  converted  to  SSH  in  order  to  be 
assimilated  by  NCOM-4DVAR.  In  this  case  a  3-yr 
mean  SSH  from  a  global  HYCOM  model  run  from 
2008  to  2011  is  used  to  estimate  the  geoid.  An  obser¬ 
vational  error  of  5  cm  for  SSHA  is  used  in  this  study;  it 
does  not  take  into  account  the  error  in  the  geoid 
estimate. 

3.  Gulf  of  Mexico  circulation  August-September 
2012 

As  mentioned  earlier,  the  GoM  is  the  region  in  which 
the  CARTHE  group  deployed  their  300  surface  drifters 
in  an  attempt  to  gather  observational  data  on  the  near¬ 
surface  circulation  of  the  Gulf.  This  is  an  interesting 
region  for  such  work,  as  the  GoM  is  a  semienclosed 
basin  that  is  dominated  by  a  strong  current  entering 
through  the  passage  between  the  Yucatan  Peninsula 
and  the  island  of  Cuba  (Yucatan  Current)  and  exiting 


through  the  Straits  of  Florida  (Florida  Current);  this 
current  is  known  as  the  Loop  Current  (LC).  The  LC 
extension  into  the  GoM  can  be  minor,  with  the  current 
moving  from  the  Yucatan  Channel  to  the  northeast 
directly  to  the  Florida  Straits  (known  as  “port  to 
port”).  At  times,  however,  the  LC  can  extend  far  into 
the  GoM,  either  to  26°N  (average  extension),  or  as  far 
north  as  28°N  (fully  extended;  Leben  2005).  When  the 
LC  is  fully  extended  a  warm-core  anticyclonic  Loop 
Current  Eddy  (LCE)  can  be  formed  by  the  completion 
of  the  circulation  around  the  northernmost  extent  of 
the  LC.  This  eddy  can  detach  and  reattach  many  times 
before  the  LCE  separates  for  good.  Once  the  LCE 
separates,  the  LC  returns  to  its  port-to-port  position 
while  the  LCE  propagates  westward  at  about  2- 
5  km  day-1.  One  possible  mechanism  for  the  detach¬ 
ment  of  the  LCE  from  the  LC  is  the  propagation  of  a 
cold-core  cyclonic  eddy  southwestward  from  the  west 
Florida  shelf,  known  as  a  Loop  Current  Frontal  Eddy 
(LCFE;  Schmitz  2005).  This  cyclonic  eddy  moves  into 
the  eastern  side  of  the  LC  and  can  penetrate  deep  to 
the  west,  eventually  resulting  in  the  detachment  of  the 
LCE  from  the  LC. 

We  can  examine  the  LC  during  the  experiment  time 
period  presented  here  using  the  V30  X  V30  mapped 
absolute  dynamic  topography  (MADT)  from  the 
Archiving,  Validation,  and  Interpretation  of  Satellite 
Oceanographic  data  (AVISO;  produced  by  Ssalto/ 
Duacs  and  distributed  by  AVISO,  with  support  from 
CNES  at  http://www.aviso.oceanobs.com/duacs/).  The 
MADT  has  been  calibrated  to  ensure  that  the  spatial 
mean  matches  that  from  the  daily  averaged  free-run 
NCOM  (initial  condition  from  global  NCOM  on  1  Au¬ 
gust,  with  no  data  assimilation  thereafter;  the  boundary 
conditions  are  listed  in  the  previous  section);  the  cali¬ 
bration  is  essentially  an  offset  to  the  MADT.  This  cali¬ 
brated  AVISO  MADT  will  be  referred  to  as  AVISO 
SSH. 

Figure  2  shows  the  GoM  SSH  for  1  August,  20  Au¬ 
gust,  10  September,  and  30  September  for  AVISO 
(Figs.  2a-d)  and  the  free-run  NCOM  (Figs.  2e-h).  The 
LC  on  1  August  seems  to  be  in  the  process  of  shedding 
a  LCE  and  returning  to  its  port-to-port  position,  with 
the  LCE  positioned  near  26°N,  88°W.  There  appears  to 
be  a  cyclonic  LCFE  located  to  the  southeast  of  the 
LCE  near  23.5°N,  85°W.  This  eddy  seems  to  be  pinch¬ 
ing  off  the  LC  from  the  LCE  as  it  extends  westward, 
and  may  be  partially  responsible  for  the  separation 
seen  in  the  AVISO  SSH.  As  time  goes  on,  from  1  Au¬ 
gust  to  30  September,  the  LCE  continues  to  remain 
unattached  to  the  LC,  with  the  cold-core  LCFE  de¬ 
caying  slowly  to  its  southeast.  The  LCE  itself  moves 
slowly  westward  and  begins  to  lose  amplitude,  going 
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Fig.  2.  Sea  surface  height  (m)  for  AVISO  (a)  1  Aug,  (b)  20  Aug,  (c)  10  Sep,  and  (d)  30  Sep;  and 
for  NCOM  free-run  (e)  1  Aug,  (f)  20  Aug,  (g)  10  Sep,  and  (h)  30  Sep. 


from  >0.5-m  height  on  1  August  to  below  0.25-m  height 
by  30  September. 

The  NCOM  free-run  solution  exhibits  a  substan¬ 
tially  different  pattern.  Figure  2e  shows  the  NCOM 
free-run  initial  condition  on  1  August.  The  LCE,  like 
what  is  seen  in  the  AVISO  SSH,  seems  to  have  de¬ 
tached  from  the  LC,  although  its  size  and  orientation 
is  slightly  different  than  the  observations.  The  LCFE 
in  the  NCOM  initial  condition  is  much  smaller  than 
what  is  seen  in  AVISO,  and  in  fact,  has  completely 
dissipated  by  20  August  (Fig.  2f).  By  10  September, 
the  LC  seems  to  extend  into  the  GoM  again,  nearly 
reattaching  with  the  LCE  (Fig.  2g).  Finally,  by 
30  September  (Fig.  2h),  the  LCE  has  now  moved 


farther  from  the  LC  while  maintaining  a  much  larger 
size  and  amplitude  than  what  is  seen  in  AVISO.  Even 
though  the  free-run  experiment  does  not  seem  to  be 
adversely  affected  by  using  the  global  NCOM  solution 
as  its  initial  condition  on  1  August  2012,  likely  due 
to  the  marginal  difference  in  model  resolutions  be¬ 
tween  the  global  model  (roughly  12  km  in  the  GoM) 
and  the  regional  model  (6km),  the  global  NCOM 
initialization  is  less  accurate  in  capturing  the  LC  and 
LCE  positions,  as  well  as  the  size,  location,  and  am¬ 
plitude  of  the  LCFE.  It  is  likely  due  to  this  that  the 
free-run  solution  diverges  quickly  from  reality  through 
the  remainder  of  the  experiment.  This  strongly  suggests 
that  an  accurate  initial  condition  is  vitally  important  to 
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Fig.  3.  (a)  NCOM  Gulf  of  Mexico  domain  with  location  of  all  assimilated  GLAD  drifter 
positions  from  1  Aug  to  30  Sep  2012  and  (b)  assimilated  GLAD  drifter  positions  during  2-6  Sep 
2012;  gray  box  indicates  the  “GLAD  region”  of  interest. 


capturing  the  evolution  of  the  LC  and  its  associated 
features  in  the  GoM. 

4.  Experiment  design  and  results 

a.  Experiment  design 

To  assess  the  impact  of  GLAD-derived  surface  velocity 
observations  on  the  model  SSH,  assimilation  experiments 
that  use  surface  velocity  observations  only  (referred  to 
here  as  VEL)  and  along-track  SSH  observations  only 
(referred  to  here  as  SSH1)  are  run,  as  well  as  the  free-run 
with  no  assimilation  (referred  to  here  as  FR);  a  final  ex¬ 
periment  that  assimilates  a  combination  of  SSH  and  ve¬ 
locity  (COM)  is  done  to  determine  their  combined  effect. 
Each  assimilation  experiment  consists  of  a  60-day  training 
period  from  1  August  to  30  September  2012;  they  all  share 
the  same  initial  condition  from  global  NCOM  on  1  August 
2012.  This  training  period  proceeds  as  a  series  of  96-h 
assimilation  windows,  where  at  the  end  of  each  window 
the  forecast  model  is  run  from  the  updated  final  condition 
to  provide  the  background  for  the  next  96-h  assimilation 
cycle.  Error  variances  and  spatial  correlation  scales  for  the 
static  covariance  applied  in  the  NCOM-4DVAR  are  set  in 
the  same  manner  as  Carrier  et  al.  (2014).  As  expected 
from  any  assimilation  experiment  attempting  to  correct  a 
significantly  wrong  model  solution  (Fig.  2),  it  takes  a  few 
assimilation  cycles  before  the  assimilated  solution  begins 
to  clearly  diverge  from  the  free  run  and  to  converge  to¬ 
ward  the  observations.  As  a  result,  the  subsequent  analysis 
shown  here  will  examine  the  solution  for  the  training 
period  from  13  August  to  30  September  (omitting  the  first 
three  assimilation  cycles  during  1-12  August  2012). 

Regarding  the  surface  velocity  observations,  the  pro¬ 
cessed  GLAD  velocities  are  available  in  15-min  in¬ 
tervals;  however,  these  data  are  only  assimilated  hourly 
due  to  the  cost  of  input/output  (I/O)  in  the  4DVAR.  The 
GLAD  data  are  thinned  at  each  hourly  time  period  by 
the  spatial  correlation  scale  (i.e.,  Rossby  radius,  roughly 


40  km  in  the  GoM)  used  in  the  4DVAR.  This  ensures 
that  no  two  assimilated  GLAD  velocity  observations  are 
within  a  correlation  scale  of  one  another.  Figure  3  shows 
the  spatial  distribution  of  the  GLAD  data  throughout 
the  60-day  experiment;  Fig.  3a  shows  the  total  obser¬ 
vation  points  through  time,  and  Fig.  3b  shows  a  sample 
distribution  of  assimilated  observations  from  one  96-h 
assimilation  window  (2-6  September  2012).  Figure  3 
shows  that  the  GLAD  drifters  are  primarily  contained 
within  the  central  GoM  for  the  duration  of  the  experi¬ 
ment;  the  western  GoM  and  the  northern  Caribbean  Sea 
are  mainly  devoid  of  observations.  The  inhomogeneous 
distribution  of  data  is  likely  to  impact  to  what  degree 
these  observations  modify  the  model  SSH.  Examining 
the  SSH  forecast  error  over  the  entire  model  domain 
would  likely  skew  the  results  and  hide  this  influence; 
therefore,  our  statistical  analysis  is  done  by  not  only 
examining  the  SSH  forecast  error  across  the  entire  do¬ 
main,  but  also  by  concentrating  on  the  area  most  densely 
covered  by  GLAD  observations,  hereafter  referred  to  as 
the  “GLAD  region,”  indicated  by  the  gray  box  in 
Fig.  3b.  The  GLAD  region  covers  the  area  of  the  LC 
where  mesoscale  eddies  are  formed  and  shed;  therefore, 
the  circulation  within  the  GLAD  region  is  the  most 
difficult  to  forecast. 

For  the  SSH1  experiment,  the  along-track  SSH  ob¬ 
servations  are  binned  and  assimilated  by  the  NCOM- 
4DVAR  every  6h  within  each  assimilation  cycle. 
Figure  4  shows  the  coverage  of  altimeter  data;  Fig.  4a 
shows  the  total  altimeter  data  assimilated  over  the  entire 
experiment  while  Fig.  4b  shows  the  coverage  during  one 
assimilation  window  (2-6  September  2012).  Figure  4 
shows  that  the  GLAD  region  is  well  covered  by  altim¬ 
eter  data  during  the  experiment,  which  should  result  in  a 
well-constrained  model  SSH  field  within  the  GLAD 
region  for  the  SSH1  experiment. 

The  96-h  forecast,  generated  from  the  4DVAR  anal¬ 
ysis  at  the  end  of  each  cycle  within  the  training  period,  is 
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compared  to  available  SSH  observations.  Along-track 
altimeter  observations  are  dense  in  the  along-track  di¬ 
rection,  but  sparse  in  space  and  time  in  the  cross-track 
direction;  therefore,  along-track  altimeter  observations 
may  not  capture  important  features  at  a  particular  in¬ 
stance  in  time.  To  avoid  this  in  the  forecast  evaluation, 
the  AVISO  SSH  has  been  selected  as  the  dataset  to 
which  the  forecasts  will  be  compared.  Using  a  composite 
SSH  map  for  forecast  validation  is  well  established  in  the 
literature  (Centurioni  et  al.  2008;  Ferron  2011;  Zavala- 
Garay  et  al.  2012;  Xu  et  al.  2013;  etc.)  and  is,  therefore, 
appropriate  for  an  examination  of  this  type. 

b.  Results  and  conclusions 

Before  examining  the  forecast  results  and  how  they 
compare  to  observations,  it  is  important  to  examine 
the  analysis  to  ensure  the  assimilation  step  is  func¬ 
tioning  properly.  With  DA  it  is  a  reasonable  require¬ 
ment  that  the  analysis  fit  the  assimilated  observations 
within  the  prescribed  observation  error.  One  can  ex¬ 
amine  this  by  calculating  a  normalized  mean  absolute 
error  that  will  be  referred  to  here  as  /fit.  The  /fit  metric 
is  calculated  as 


'fit 


i  M 

=-  y 


I  y  -H  xa  I 

m  m  1 


(10) 


where  xa  is  the  model  analysis,  ym  is  the  rath  observation, 
crm  is  the  observation  error,  H  is  an  operator  that  maps 


the  state  variables  to  the  observation  space,  and  M  is  the 
total  number  of  observations.  From  (10)  it  can  be  shown 
that  the  analysis  is  fitting  the  observations  within  the 
prescribed  observation  error  if  the  value  of  /fit  is  at  or 
less  than  1.0.  Figure  5  shows  the  /fit  values  for  the  SSH 
for  the  first  guess  (i.e.,  the  background  96-h  forecast,  FG; 
gray)  and  the  analysis  (AN;  black)  for  the  SSH1  experi¬ 
ment  from  13  August  to  30  September  2012  (as  compared 
against  assimilated  along-track  SSH  observations).  Each 
/fit  value  is  computed  by  comparing  the  first  guess  and 
analysis,  valid  at  0000  UTC  daily,  to  available  SSH  ob¬ 
servations  that  fall  ±3  h  around  this  0000  UTC  time  level. 
Figure  5  indicates  that  the  NCOM-4DVAR  is  generally 
fitting  the  along-track  SSH  observations  within  error 
limits  for  the  entire  experiment  time  frame  (the  value  of 
1.0  is  indicated  by  the  dashed  line). 

This  check  can  be  performed  for  the  VEL  experiment 
as  well.  Figure  6  shows  the  /fit  value  for  the  surface  ve¬ 
locity  field  for  the  first  guess  (gray)  and  the  analysis 
(black)  for  the  VEL  experiment  from  13  August  to 
30  September  2012  (cf.  assimilated  surface  velocity  ob¬ 
servations).  As  with  the  SSH1  experiment,  the  VEL 
analysis  is  generally  fitting  the  observations  within  the 
prescribed  error.  There  is  one  point  in  late  August 
where  this  is  not  the  case,  as  /fit  values  during  this  time 
are  near  2.0.  This  coincides  with  the  arrival  of  Hurricane 
Isaac  over  the  northern  GoM,  with  the  storm  track  di¬ 
rectly  through  the  GLAD  region.  Figure  7  shows  the 
NOGAPS  surface  wind  stress  magnitude  (interpolated 


Fig.  5.  The  /fit  metric  plot  for  SSH  from  the  SSH1  experiment. 
First  guess  (FG)  field  is  in  gray  and  analysis  (AN)  is  in  black. 
Values  computed  by  comparison  to  assimilated  SSH  observations 
from  13  Aug  to  30  Sep  2012. 


Fig.  6.  The  /fit  metric  plot  for  velocity  from  the  VEL  experiment 
as  compared  to  assimilated  GLAD  drifter-derived  surface  veloci¬ 
ties  (13  Aug-30  Sep  2012). 
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Fig.  7.  NOGAPS  surface  wind  stress  (interpolated  to  the  NCOM  domain)  for  (a)  28  Aug 
and  (b)  29  Aug  2012.  Wind  stress  direction  is  indicated  by  arrows  and  magnitude  is  indicated 
by  color  contours  (in  N  m-2).  GLAD  region  is  indicated  by  a  gray  box. 


to  the  NCOM  grid)  for  28  and  29  August  (Figs.  7a  and 
7b,  respectively)  with  wind  direction  indicated  by  curved 
arrows;  Fig.  7  also  shows  the  GLAD  region  as  a  gray  box 
outline.  By  28  August,  the  NOGAPS  simulation  has 
Isaac  entering  the  GLAD  region  and  making  landfall 
later  on  29  August.  Hurricanes  produce  very  strong 
surface  wind  forcing  that  alters  surface  ocean  currents 
more  toward  ageostrophy.  In  addition  to  this,  the  sur¬ 
face  wind  forcing  for  NCOM  is  from  the  0.5°  NOGAPS 
model,  a  resolution  that  does  not  allow  the  NOGAPS 
model  to  fully  capture  the  intensity  of  the  surface  wind. 
For  example,  on  28  August  (at  0000  UTC),  NOGAPS- 
simulated  peak  storm  winds  were  19.6 ms-1,  whereas 
the  National  Hurricane  Center  listed  Isaac  with  winds 
near  31  ms-1.  Also,  at  this  resolution  NOGAPS  is  not 
capable  of  resolving  small-scale  features  in  those  winds 


that  likely  have  a  significant  impact  on  the  surface 
drifters  from  GLAD.  Thus,  the  background  ocean 
state  used  in  the  assimilation  does  not  capture  the 
realistic  ocean  surface  currents  during  this  extreme 
weather  event.  Further,  the  prescribed  model  and 
observation  errors  used  here  in  NCOM-4DVAR  do 
not  change  with  time.  These  cumulative  factors  are 
likely  responsible  for  a  poor  fit  to  the  GLAD  obser¬ 
vations  in  the  analysis  during  this  time  period.  Once 
the  storm  passed,  however,  the  analysis  fit  to  the  sur¬ 
face  velocity  observations  improves  significantly  and 
rapidly,  as  the  /fit  values  fall  once  again  to  and  below 
1.0  for  the  remainder  of  the  training  period.  Also,  the 
first-guess  field  exhibits  a  downward  trend  in  error, 
beginning  with  a  /fit  value  of  3.0  on  13  August  and 
ending  with  a  value  just  above  2.0  by  30  September. 


Fig.  8.  RMS  error  over  (a)  entire  GoM  domain  and  (c)  GLAD  region  only;  correlation  over 
(b)  entire  GoM  domain  and  (d)  GLAD  region  only.  Statistics  shown  for  FR  (thin  black),  SSH1 
(thick  black),  and  VEL  (gray)  experiment  96-h  SSH  forecast  solutions  (cf.  available  AVISO 
SSH  observations).  SSH  observation  error  indicated  in  (a)  and  (c)  by  dashed  black  line. 
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Fig.  9.  Sea  surface  height  (m)  for  AVISO  (a)  1  Aug,  (b)  20  Aug,  (c)  10  Sep,  and  (d)  30  Sep;  for  SSH1  experiment 
(e)  1  Aug,  (f)  20  Aug,  (g)  10  Sep,  and  (h)  30  Sep;  and  for  VEL  experiment  (i)  1  Aug,  (j)  20  Aug,  (k)  10  Sep,  and 
(1)30  Sep. 


This  indicates  that  the  background  state  is  improving 
with  time  as  the  model  continues  to  train  toward  ob¬ 
servations.  The  fact  that  this  downward  trend  is  not 
seen  with  the  SSH  assimilation  (Fig.  5)  is  not  sur¬ 
prising,  as  the  placement  of  along-track  SSH  obser¬ 
vations  from  one  time  to  the  next  can  be  very 
different,  whereas  the  GLAD  observations  are  pri¬ 
marily  confined  to  the  central  Gulf  throughout  the 
experiment. 

An  examination  can  be  made  of  the  model  fit  to 
AVISO  during  the  training  period  from  the  FR, 
SSH1,  and  VEL  experiments  in  order  to  ascertain  if 
the  model  solution  from  the  assimilation  experiments 
is  improved  over  the  free  run.  To  do  this  the  96-h 
forecast  generated  from  each  analysis  during  the 
training  period  is  compared  to  the  AVISO  SSH.  The 


root-mean-square  (RMS)  error  and  correlation  are 
given  as 


/ 1  M 

RMS  =  W—  X  (y  ~H  xrf  and 

l  /  /x  /l  m  m  ; 


(li) 


r  =  - 


l(Hmxf-X)(ym-Y) 

m= 1 


M 


X  (Hmxf  -  x  f 

m= 1 


M 

^(ym-yy 

m= 1 


:  ,  (12) 


where  is  the  model  forecast,  X  is  the  model  mean  at 
the  observation  locations  (as  mapped  by  H ),  and  Y  is  the 
mean  of  the  observations;  and  r  is  known  as  the  Pearson 
product-moment  correlation  coefficient.  To  statistically 


March  2016 


CARRIER  ET  AL. 


1061 


a  b 


c 


d 


0.2 


Fig.  10.  (a)  RMS  error  and  (b)  correlation  over  entire  GoM  domain;  (c)  RMS  error  and 
(d)  correlation  within  the  GLAD  region.  Statistics  shown  for  SSH1  (solid  black),  VEL  (solid 
gray),  SSH1  persistence  (dashed  black),  and  VEL  persistence  (dashed  gray)  SSH  30-day 
forecast  trajectory  experiments  (cf.  available  AVISO  SSH  observations). 


compare  to  AVISO,  the  NCOM  SSH  fields  are  a  daily 
average  of  the  last  day  of  each  96-h  forecast.  Figure  8 
shows  the  RMS  error  (Figs.  8a  and  8c)  and  correlation 
(Figs.  8b  and  8d)  for  the  FR  (thin  black),  SSH1  (thick 
black),  and  VEL  (gray)  experiments  over  the  whole 
model  domain  (top  panels)  and  just  within  the  GLAD 
region  (bottom  panels);  Figs.  8a  and  8c  also  show  the 
prescribed  SSH  observation  error  (dashed  black)  for 
comparison.  There  is  some  improvement  seen  in  the 
SSH1  and  VEL  experiments  over  the  FR  when  exam¬ 
ining  the  RMS  error  over  the  entire  domain  (Fig.  8a),  as 
the  FR  error  tops  out  at  about  0.15  m,  while  the  SSH1 
and  VEL  error  stay  just  at  or  slightly  above  0.1  m  by  the 
end  of  the  training  period.  The  SSH1  seems  to  out¬ 
perform  the  VEL  experiment  in  early  September, 
however,  the  error  in  both  SSH1  and  VEL  when  com¬ 
pared  to  AVISO  becomes  nearly  identical  by  the  end  of 
the  training  period.  Likewise,  the  correlation  to  AVISO 
over  the  entire  domain  is  also  quite  comparable 
(Fig.  8b),  with  the  SSH1  experiment  showing  slightly 
better  correlation  with  AVISO  than  the  VEL  experi¬ 
ment.  When  focusing  in  on  the  GLAD  region,  however, 
there  is  a  stark  improvement  in  the  SSH1  and  VEL 
forecasts  over  the  FR  (Fig.  8c).  The  FR  error  reaches  as 
high  as  0.21  m,  while  the  error  in  SSH1  and  VEL  is  near 
0.15  m  around  the  time  of  Hurricane  Isaac’s  passing  (28- 
29  August)  then  drops  to  about  0.1  m  by  the  end  of  the 
training  period.  Within  the  GLAD  region,  the  VEL 
model  forecast  fits  the  AVISO  SSH  just  as  well  as  the 
SSH1  model  forecast;  this  indicates  the  usefulness  of 
the  velocity  observations  when  attempting  to  constrain 


the  model  SSH  field.  The  model  SSH  forecast  correla¬ 
tion  to  AVISO  data  within  the  GLAD  region  is  also 
telling:  the  VEL  experiment  shows  slightly  higher  cor¬ 
relation  than  the  SSH1  experiment,  especially  through 
the  month  of  September.  This  is  likely  due  to  the  fact 
that  the  VEL  experiment  has  many  more  observations, 
of  an  in  situ  nature,  within  the  GLAD  region  than  does 
the  SSH1  experiment. 

A  qualitative  comparison  between  the  FR,  SSH1,  and 
VEL  SSH  forecasts  can  be  made  by  comparing  the 
model  solutions  to  each  other  and  to  AVISO.  Figure  9 
shows  a  comparison  of  the  SSH  from  AVISO  on  1  Au¬ 
gust  (Fig.  9a),  20  August  (Fig.  9b),  10  September 
(Fig.  9c),  and  30  September  (Fig.  9d),  to  the  corre¬ 
sponding  SSH1  experiments  (Figs.  9e-h)  and  VEL  ex¬ 
periments  (Figs.  9i— 1).  As  was  shown  in  Fig.  2,  the  initial 
condition  obtained  from  global  NCOM  (Figs.  9e,i) 
shows  some  significant  differences  when  compared  to 
AVISO  (Fig.  9a).  It  was  shown  how  the  free-run  NCOM 
never  matches  AVISO  throughout  the  experiment  time 
frame,  but  with  the  assimilation  of  SSH  or  velocity  ob¬ 
servations,  the  results  are  quite  different.  Figure  9f 
shows  the  SSH1  forecast  on  20  August  and  already  the 
SSH  field  is  beginning  to  match  quite  well  with  AVISO. 
The  warm-core  LCE  eddy  (26°N,  88°W)  is  elongated  in 
the  northwest-southeast  direction,  very  similarly  to 
AVISO,  and  the  cold-core  LCFE  to  the  southeast  of  the 
LCE  seen  in  AVISO  is  now  apparent  in  the  SSH1  ex¬ 
periment  (unlike  the  FR).  The  SSH1  experiment,  how¬ 
ever,  does  not  exhibit  the  “bow”  shape  in  the  LCE  on 
20  August,  likely  due  to  a  weaker  secondary  cold-core 
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Fig.  11.  Sea  surface  height  (m)  for  AVISO  (a)  1  Oct,  (b)  15  Oct,  and  (c)  30  Oct;  for  SSH1  experiment  (d)  1  Oct,  (e)  15  Oct,  and  (f)  30  Oct; 

and  for  VEL  experiment  (g)  1  Oct,  (h)  15  Oct,  and  (i)  30  Oct. 


eddy  to  the  northeast  of  the  LCE  (seen  in  AVISO, 
Fig.  9b  near  27°N  and  87°W).  The  VEL  experiment,  on 
the  other  hand,  does  seem  to  capture  the  shape  of  the 
LCE  better  on  20  August  (Fig.  9j),  albeit  with  greater 
amplitude  at  the  eddy  center.  The  VEL  experiment  does 
seem  to  capture  both  cold-core  eddies,  however,  the 
amplitudes  are  lower  than  what  is  seen  in  AVISO.  The 
SSH1  experiment  seems  to  match  quite  well  with 
AVISO  on  10  September  (Figs.  9c, g),  with  a  dual  max¬ 
imum  seen  in  the  LCE  and  the  cold-core  LCFE  to  its 
southeast.  However,  the  SSH1  experiment  exhibits  an 
LCFE  (24°N,  86°W)  that  is  deeper  than  what  is  seen  in 
AVISO  and  it  seems  to  connect  this  eddy  with  the  sec¬ 
ondary  cold-core  eddy  (seen  in  Fig.  9c  to  the  northeast  of 
the  primary  LCE).  Here  the  SSH1  experiment,  unlike 
the  FR  (Fig.  2f),  does  not  extend  the  LC  farther  into  the 
GoM,  and  the  LCE  is  not  close  to  reattachment.  This  is 
likely  due  to  the  stronger  LCFE  in  the  SSH1  experiment 
than  what  the  FR  exhibits.  The  VEL  experiment  on 
10  September  (Fig.  9k)  shows  a  more  smoothed  LCE 


than  what  is  shown  in  AVISO  or  the  SSH1  experiment 
at  this  time,  though  it  is  much  closer  to  AVISO  and 
SSH1  than  the  FR  (Fig.  2g).  The  VEL  experiment  seems 
to  have  the  dual  maximum  in  the  LCE,  as  in  AVISO, 
though  with  a  higher  core  amplitude.  Also,  the  VEL 
experiment  does  exhibit  the  cold-core  LCFE,  near  24°N, 
86°W,  though  it  is  weaker  than  what  is  shown  in  AVISO 
and  SSH1.  By  the  end  of  the  training  period,  30  Sep¬ 
tember,  the  LCE  in  AVISO  seems  to  be  decaying,  as  the 
core  amplitude  has  dropped  significantly  since  1  August 
(Fig.  9d).  The  SSH1  experiment  (Fig.  9h)  generally  ex¬ 
hibits  the  same  structure  of  the  LCE  and  LCFE  as  seen 
in  AVISO,  but  again  with  slightly  larger  amplitudes  in 
both  eddies.  The  VEL  experiment  (Fig.  91),  on  the  other 
hand,  not  only  matches  the  general  structure  but  is 
closer  in  amplitude  to  AVISO  as  well.  It  is  important  to 
note  that  the  assimilation  of  either  SSH  or  velocity  ob¬ 
servations  seems  to  constrain  the  SSH  field  well  (when 
compared  to  AVISO)  and  also  helps  to  capture  the  ex¬ 
istence  of  the  cold-core  LCFE  to  the  southeast  of  the 
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Fig.  12.  RMS  error  over  (a)  entire  GoM  domain  and  (c)  GLAD  region  only;  correlation  over 
(b)  entire  GoM  domain  and  (d)  GLAD  region  only.  Statistics  shown  for  SSH1  (thick  black), 
VEL  (gray),  and  COM  (dashed  gray)  experiment  96-h  SSH  forecast  solutions  (cf.  available 
AVISO  SSH  observations).  SSH  observation  error  indicated  in  (a)  and  (c)  by  dashed  black  line. 


LCE,  which  was  generally  missed  by  the  NCOM  FR 
experiment. 

It  has  been  established  that  the  assimilation  of  the 
GLAD -derived  surface  velocity  observations  is  useful  in 
constraining  the  model  SSH  field,  nearly  as  well  as  the 
assimilation  of  along-track  SSH  itself.  Given  this,  it 
would  be  useful  to  examine  the  forecast  SSH  error  for 
each  experiment  for  a  longer  time  period  than  the  96  h 
that  has  been  shown.  Examining  the  forecast  error  out  to 
30  days,  for  instance,  should  shed  light  on  the  extent  to 
which  the  model  SSH  is  constrained. 

For  this  examination,  30-day  forecasts  are  run  from 
the  final  analysis  state  on  30  September  2012  produced 
by  the  training  period  for  the  SSH1  and  VEL  experi¬ 
ments  (using  the  same  surface  and  lateral  boundary 
forcing  sources  as  described  in  section  2).  The  model 
forecast  SSH  is  compared,  as  before,  to  the  available 
AVISO  SSH  observations.  Figure  10  shows  the  RMS 
error  (Fig.  10a)  and  correlation  (Fig.  10b)  for  the  SSH1 
(thick  black)  and  VEL  (gray)  experiments  for  the  en¬ 
tire  model  domain;  Figs.  10c  and  lOd  show  the  RMS 
error  and  correlation,  respectively,  for  the  GLAD  re¬ 
gion.  The  persistence  forecast  (i.e.,  nonevolving  model 
state  from  30  September)  is  also  shown  for  SSH1 
(dashed  black;  SSHl-p)  and  VEL  (dashed  gray;  VEL-p) 
experiments.  Examining  the  results  in  the  full  do¬ 
main  (Figs.  10a, b),  it  is  clear  that  the  forecast  error  for 
both  SSH1  and  VEL  grows  slowly  through  time,  from 
0.1  m  on  30  September  to  just  around  0.12  by  30  October 
(correlation  dropping  from  0.9  to  near  0.8).  The  persis¬ 
tence  forecast  for  SSH1  grows  in  error  (from  0.1  to 


0.15  m)  more  steeply  than  the  SSH1  forecast  (and  cor¬ 
relation  drops  sharply,  from  0.9  to  below  0.8),  indicating 
some  measure  of  skill  for  the  SSH1  model  forecast.  The 
same  can  be  said  for  the  VEL  experiment,  though  its 
persistence  forecast  error  grows  less  sharply.  The  results 
in  the  GLAD  domain  (Figs.  10c, d)  are  roughly  the  same, 
with  both  the  SSH1  and  VEL  experiments  exhibiting 
forecast  skill  over  persistence.  The  FR  experiment  (not 
shown)  exhibits  much  higher  error  in  both  the  forecast 
and  persistence  than  either  SSH1  or  VEL  through  all  of 
October. 

The  30-day  forecasts  can  be  examined  qualitatively  as 
well  by  comparing  the  SSH  fields  to  each  other  and  to 
AVISO.  Figure  11  shows  the  SSH  field  from  AVISO  on 
1  October  (Fig.  11a),  15  October  (Fig.  lib),  and  30  Oc¬ 
tober  (Fig.  11c);  also  shown  are  the  forecast  SSH  field 
from  the  SSH1  experiment  on  1  October  (Fig.  lid), 
15  October  (Fig.  lie),  and  30  October  (Fig.  Ilf),  and 
from  the  VEL  experiment  on  1  October  (Fig.  llg), 
15  October  (Fig.  llh),  and  30  October  (Fig.  lli).  Over 
the  month  of  October  one  can  see  that  the  primary  LCE 
continues  to  decay  as  it  moves  westward  in  the  AVISO 
maps  (Figs,  lla-c);  also,  by  30  October,  the  AVISO 
observations  indicate  that  the  LC  is  beginning  to  slightly 
extend  back  into  the  GoM.  The  SSH1  forecast  SSH  field 
on  1  October  (Fig.  lid)  seems  close  to  that  of  AVISO; 
the  primary  LCE  is  near  25°N,  89°W  with  an  accompa¬ 
nying  cold-core  eddy  to  its  southeast  at  24°N,  86°W. 
Though,  it  should  be  noted  that  the  SSH1  forecast  LC 
and  cold-core  eddies  both  have  higher  amplitudes  than 
what  is  suggested  by  AVISO.  The  VEL  forecast  SSH 
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Fig.  13.  Sea  surface  height  (m)  for  AVISO  (a)  1  Aug,  (b)  20  Aug,  (c)  10  Sep,  and  (d)  30  Sep;  and 
for  COM  experiment  (e)  1  Aug,  (f)  20  Aug,  (g)  10  Sep,  and  (h)  30  Sep. 


field  on  1  October  (Fig.  llg)  appears  to  be  closer  to 
AVISO  than  the  SSH1  experiment.  The  VEL  experi¬ 
ment  LC  and  cold-core  eddies  are  positioned  in  similar 
locations  as  those  found  in  the  SSH  experiment;  how¬ 
ever,  the  amplitudes  of  these  features  are  slightly  lower 
in  the  VEL  experiment  when  compared  to  the  SSH 
experiment  (especially  for  the  cold-core  eddy).  This 
pattern  continues  to  15  October  (Figs.  llb,e,h)  as 
the  SSH1  forecast  SSH  continues  to  exhibit  higher- 
amplitude  LC  and  cold-core  eddies  than  what  is  seen 
in  AVISO;  meanwhile,  the  VEL  SSH  field  shows  a 
quickly  decaying  cold-core  eddy  to  the  southeast  of  the 
LCE,  which  still  exhibits  higher  amplitude  than  AVISO. 
Finally,  on  30  October  (Figs.  llc,f,i),  AVISO  continues 


to  show  the  LCE  in  a  more  decayed  state  than  either  the 
SSH1  or  VEL  model  results.  The  SSH1  experiment 
continues  to  exhibit  a  relatively  strong  cold-core  eddy  to 
the  east  of  the  LCE,  which  seems  to  be  suppressing 
any  intrusion  of  the  LC  unlike  what  is  seen  in  the 
AVISO  observations.  The  VEL  experiment,  while 
still  exhibiting  a  higher- amplitude  LCE,  shows  that  the 
cold-core  eddy  is  nearly  decayed  and,  as  a  possible 
consequence,  the  LC  seems  to  be  intruding  somewhat 
northward  into  the  southern  GoM  similar  to  AVISO. 
The  SSH  fields  shown  here  support  the  statistical  find¬ 
ings  shown  in  Fig.  10  that  suggest  that  the  VEL  forecast, 
while  exhibiting  slightly  higher  error  than  the  SSH1 
forecast  early  in  October,  seems  to  be  performing  as  well 
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Fig.  14.  NCOM  surface  velocity  forecast  RMS  error  for  the  SSH1  (thick  black)  and  COM 
(dashed  gray)  experiments  (computed  using  derived  Eulerian  velocities  from  GLAD  drifter 
observations).  Markers  indicate  the  start/end  of  each  96-h  forecast  cycle. 


as,  if  not  better  than,  the  SSH1  forecast.  This  suggests 
that  the  forecast  error  in  the  VEL  and  SSH1  experi¬ 
ments  are  roughly  identical,  indicating  that  the  assimi¬ 
lation  of  velocity  measurements  has  helped  to  constrain 
the  model  SSH  in  a  comparable  fashion  with  the  SSH1 
experiment. 

Finally,  it  should  be  determined  how  the  analysis  and 
resulting  forecast  performs  when  both  along-track  SSH 
and  GLAD -derived  velocity  observations  are  assimi¬ 
lated.  This  gives  a  sense  of  the  added  value  from  the 
inclusion  of  velocity  observations  with  the  standard 
set  of  temperature,  salinity,  and  SSH  observations. 
Figure  12  shows  the  RMS  (Figs.  12a, c)  and  the  correla¬ 
tion  of  each  of  the  96-h  forecasts  during  the  training 
period  (Figs.  12b, d)  for  the  SSH1  (thick  black)  and  VEL 
(gray)  experiments,  as  well  as  the  combined  assimilation 
(COM,  dashed  gray)  for  the  entire  Gulf  region  (top 
panels)  and  the  GLAD  region  (bottom  panels)  as 
compared  to  AVISO  SSH;  Figs.  12a  and  12c  also  show 
the  prescribed  SSH  observation  error  (dashed  black) 
for  comparison.  Generally,  the  COM  SSH  forecast 
performs  closely  with  the  SSH1  experiment  throughout 
the  Gulf  (Figs.  12a  and  12b),  but  outperforms  both  the 
SSH1  and  VEL  experiments  within  the  GLAD  region 
(Figs.  12c  and  12d).  This  indicates  that  the  observation 
types  (i.e.,  GLAD  velocities  and  along-track  SSH)  are 


complementary  to  each  other.  This  is  likely  because  the 
GLAD  velocities  act  to  “fill  in”  the  gaps  in  the  satellite 
coverage  from  the  altimeters  within  the  GLAD  region. 
This  allows  for  a  better  analysis  of  both  the  SSH  and 
surface  velocity  fields  and,  therefore,  a  superior  forecast 
as  compared  to  assimilating  either  observation  type 
alone.  Figure  13  shows  the  qualitative  comparison  of  the 
SSH  field  on  1  August,  20  August,  10  September,  and 
30  September  for  AVISO  (Figs.  13a-d)  and  the  COM 
experiment  (Figs.  13e-h).  The  COM  experiment  seems 
to  capture  the  primary  features  seen  in  AVISO  better 
than  the  SSH1  or  VEL  experiments  alone;  for  instance, 
on  20  August  (Figs.  13b  and  13f),  the  COM  experiment 
not  only  captures  the  LCFE  better  than  the  VEL  ex¬ 
periment,  but  it  also  captures  the  elongated  “bow” 
shape  to  the  LCE  that  was  not  captured  by  the  SSH1 
experiment.  And,  on  10  September  (Figs.  13c  and  13g), 
the  COM  experiment  simulates  the  cold-core  LCFE 
slightly  better  than  the  SSH1  experiment,  but  signifi¬ 
cantly  better  than  the  VEL  experiment,  while  still  sim¬ 
ulating  the  dual  maximum  in  the  LCE.  Finally,  on 
30  September  (Figs.  13d  and  13h),  the  COM  experiment 
seems  to  capture  the  lower  amplitude  and  orientation  of 
the  LCE,  while  also  matching  the  LCFE  in  AVISO 
better  than  in  the  SSH1  experiment.  It  appears  that  the 
assimilation  of  GLAD -derived  velocities  has  helped  in 


92°W  90°W  88°W  86°W  84°W  82°W 


Fig.  15.  The  96-h  forecast  SSH  (color  contours)  and  model  drifter  trajectories  (purple) 
compared  with  observed  GLAD  drifter  trajectories  (green)  for  (a)  SSH1  and  (b)  COM  ex¬ 
periments.  Forecast  is  valid  during  21-25  Aug  2012  (SSH  field  shown  is  model  average  over 
21-25  Aug). 
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the  COM  experiment  to  reduce  the  amplitude  of  the 
LCFE  from  what  is  seen  in  the  SSH1  experiment. 
Overall,  the  combination  of  along-track  SSH  and  ve¬ 
locity  observations  provides  the  best  forecast  compared 
to  AVISO  SSH  than  either  the  SSH1  or  VEL  experi¬ 
ments  alone. 

A  similar  examination  can  be  made,  but  this  time 
focusing  on  the  surface  velocity  rather  than  the  SSH 
forecast.  This  can  be  done  by  comparing  the  SSH1  and 
COM  surface  velocity  forecasts  against  the  GLAD- 
derived  surface  velocity  observations.  It  is  true  that  the 
GLAD  observations  are  assimilated  in  COM,  but  as 
noted  in  Carrier  et  al.  (2014),  comparing  the  96-h 
forecast  trajectories  to  observations  can  be  consid¬ 
ered  an  independent  examination  as  these  velocities 
have  not  yet  been  assimilated.  Figure  14  shows  the 
RMS  error  for  the  SSH1  (thick  black)  and  COM 
(dashed  gray)  surface  velocity  forecasts  (line  markers 
indicate  beginning/end  of  each  96-h  forecast  trajec¬ 
tory).  The  COM  experiment  demonstrates  superior  fit 
to  the  GLAD  velocity  observations  when  compared  to 
the  SSH1  experiment,  indicating  that  the  inclusion  of 
velocity  observations  in  the  assimilation  helps  to  fur¬ 
ther  constrain  the  model  velocity  field  better  than  SSH 
observations  alone.  This  apparent  difference  in  the 
quality  of  the  forecast  velocities  can  be  seen  more 
clearly  when  examining  model  drifter  trajectories.  In 
this  case,  model  drifters  are  initialized  at  the  same  lo¬ 
cations  as  actual  observed  GLAD  drifters.  The  model 
drifters  are  then  advected  with  the  model  forecast  ve¬ 
locity  fields  through  the  96-h  forecast.  The  trajectories 
over  this  time  can  be  compared  to  the  actual  observed 
trajectories  from  the  GLAD  drifters,  as  was  done  in 
Muscarella  et  al.  (2015).  Muscarella  et  al.  (2015)  point 
out  that  small  discrepancies  in  surface  velocity  can  re¬ 
sult  in  large  differences  in  drifter  trajectories.  Figure  15 
shows  the  model  (purple)  versus  GLAD  (green)  drifter 
trajectories  over  one  96-h  forecast  time  period  (21- 
25  August,  2012)  from  the  SSH1  (Fig.  15a)  and  the 
COM  experiments  (Fig.  15b);  initial  drifter  positions 
are  indicated  by  black  crosses.  The  trajectories  are 
overlaid  on  the  average  NCOM  SSH  field  (color  con¬ 
tours)  from  21  to  25  August.  Figure  15  clearly  demon¬ 
strates  that  the  forecast  drifter  trajectories  from  the 
COM  experiment  match  the  GLAD  trajectories  far 
more  closely  than  from  the  SSH1  experiment.  This  is 
likely  due  to  the  effect  that  the  additional  velocity 
observations  have  on  the  forecast  SSH  structure,  as  the 
COM  SSH  field  exhibits  a  slight  “L  shaped”  pattern  in 
the  primary  LCE  (26°N,  88°W)  that  is  not  seen  in  the 
SSH1  experiment.  Here,  the  velocity  observations  help 
to  fill  in  the  gaps  between  the  altimeter  tracks  and 
provide  a  more  complete  estimate  of  the  underlying 


mesoscale  structure  to  the  analysis,  which  results  in  the 
improved  forecast. 

5.  Summary 

An  examination  of  the  impact  of  GLAD -derived 
surface  velocity  observations  on  the  model  sea  surface 
height  in  the  GoM  is  presented.  These  results  are  com¬ 
pared  to  a  free-run  (i.e.,  no  assimilation)  model,  to  a 
model  that  assimilates  the  along-track  SSH  observations 
directly,  and  to  AVISO  SSH  maps.  The  results  pre¬ 
sented  here  indicate  that  while  the  free-run  model 
exhibits  error  growth  that  strongly  impacts  eddy  po¬ 
sitioning  and  amplitude  by  the  end  of  the  experiment 
time  frame,  the  assimilation  of  GLAD -derived  surface 
velocity  observations  can  help  to  constrain  model  SSH 
features  in  a  manner  consistent  with  the  assimilation  of 
along-track  SSH  observations.  It  is  further  shown  that  a 
month-long  forecast  generated  from  an  analysis  that 
utilized  the  cycled  assimilation  of  GLAD-derived  sur¬ 
face  velocities  exhibits  very  similar  forecast  error  growth 
when  compared  to  a  30-day  forecast  generated  from  an 
analysis  that  used  along-track  SSH  observations.  Fi¬ 
nally,  it  is  demonstrated  that  the  assimilation  of  both 
along-track  SSH  and  velocities  provides  the  best  fore¬ 
cast  compared  to  AVISO  SSH,  suggesting  that  the 
velocity  observations  are  complementary  to  the  along- 
track  SSH  and  help  to  “fill  in”  the  coverage  gaps  of 
the  altimeter  observations,  which  is  especially  clear 
when  examining  the  model-observed  drifter  trajectory 
comparisons. 

It  should  be  noted  that  obtaining  corrections  to  the  sea 
surface  height  from  velocity  measurements  is  only  at¬ 
tainable  when  using  a  sufficiently  sophisticated  DA 
system,  such  as  the  4DVAR.  Less-sophisticated  assimi¬ 
lation  methods,  such  as  the  3DVAR,  may  have  trouble 
accounting  for  the  time-varying  nature  of  the  velocity 
observations,  as  well  as  the  dynamical  balance  re¬ 
lationships  that  are  required.  Therefore,  in  order  to 
obtain  the  maximum  amount  of  useful  information  from 
velocity  observations,  one  must  account  for  the  tempo¬ 
ral  nature  of  the  observations  as  well  as  account  for  the 
cross  covariances  to  each  of  the  other  model  variables. 

It  has  now  been  demonstrated  that  surface  velocity 
observations  not  only  help  to  constrain  the  model  ve¬ 
locity  field,  but  can  also  be  used  to  correct  the  model 
SSH — a  very  important  necessity,  especially  given  the 
sparse  nature  of  altimeter  observations.  This  study 
suggests  that  targeted  drifter  observations  can  be  used  to 
help  constrain  the  model  SSH  with  applications  in,  not 
only  large-scale  forecasting,  but  also  in  regional  appli¬ 
cations  (e.g.,  capturing  the  location  of  the  loop  current 
edge  for  oil  rig  operations). 
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